MODULE SRFVEGEVOL_MOD
CONTAINS
SUBROUTINE SRFVEGEVOL(KIDIA,KFDIA,KLON,KLEVS,KSTEP, &
 & PTSPHY, PLAT,PMU0, &
 & PCVT,PCVL,PCVH, KVEG,ITVL,ITVH, &
 & PTSKM1M,PTSAM1M, &
 & PANFMVT, PANDAYVT, &
 & PLAIVT,PLAIL,PLAIH, &
 & PLAILC,PLAIHC, &
 & PRESPBSTR,PRESPBSTR2,PBIOMASS_LAST,PBIOMASSTR_LAST,PBIOMASSTR2_LAST, &
 & PBLOSSVT, PBGAINVT, & 
 & PLAIE1, PBSTRE1, PBSTR2E1, &
 & PLAI, PBIOM, PBLOSS, PBGAIN, PBIOMSTR, PBIOMSTR2, &
 & PDHBIOS,PDHVEGS, &
 & PSSDP2, &
 & YDCST,YDVEG,YDAGS)


! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**  *SRFVEG_EVOL*  CALLED BY *SURFTSTP_CTL*

!     PURPOSE
!     -------
!     performs the time evolution of vegetation parameters
!     at solar midnight in the case of interactive vegetation 

!     PARAMETER     DESCRIPTION                                   UNITS
!     ---------     -----------                                   -----
!     INPUT PARAMETERS (INTEGER):

!    *KIDIA*        START POINT
!    *KFDIA*        END POINT
!    *KLON*         NUMBER OF GRID POINTS PER PACKET
!    *KLEVS*
!    *KSTEP*        CURRENT TIME STEP INDEX

!     INPUT PARAMETERS (REAL):

!    *PTSPHY*       TIME STEP FOR THE PHYSICS
!    *PLAT*         LATITUDE
!    *PCVT*         VEGETATION TYPE FRACTION                      (0-1)
!    *PTSKM1M*      SKIN TEMPERATURE                              K
!    *PTSAM1M*      SOIL TEMPERATURE                              K
!    *PANFMVT*      MAXIMUM LEAF ASSIMILATION              KG_CO2/KG_AIR M/S    
!                   positive downwards
!    *PANDAYVT*     DAILY NET CO2 ASSIMILATION OVER CANOPY    KG_CO2/M2
!                   positive downwards

!     UPDATED PARAMETERS (REAL):

!    *PLAIVT*       LEAF AREA INDEX                           (M2/M2)
!    *PRESPBSTR*    RESPIRATION OF ABOVE GROUND STRUCTURAL BIOMASS    KG_CO2/M2
!    *PRESPBSTR2*   RESPIRATION OF BELOW GROUND STRUCTURAL BIOMASS    KG_CO2/M2
!    *PBIOMASS_LAST* (ACTIVE) LEAF BIOMASS OF PREVIOUS DAY               KG/M2
!     NB: value only after nitro_decline, not after laigain!!!!
!    *PBIOMASSTR_LAST* ABOVE GROUND STRUCTURAL BIOMASS OF PREVIOUS DAY   KG/M2
!    *PBIOMASSTR2_LAST* BELOW GROUND STRUCTURAL BIOMASS OF PREVIOUS DAY  KG/M2
!    *PBLOSSVT*     ACTIVE BIOMASS LOSS                                  KG/M2
!    *PBGAINVT*     ACTIVE BIOMASS GAIN                                  KG/M2
!    *PLAIE1*       TENDENCY OF LAI                                  (M2/M2)/S
!    *PBSTRE1*      TENDENCY OF ABOVE GROUND STRUCTURAL BIOMASS          KG/M2
!    *PBSTR2E1*     TENDENCY OF BELOW GROUND STRUCTURAL BIOMASS          KG/M2

!     OUTPUT PARAMETERS (REAL):

!    *PLAI*         GRID AVERAGED LAI                         (M2/M2)
!    *PBIOM*        GRID AVERAGED BIOMASS                     (KG/M2)
!    *PBLOSS*       GRID AVERAGED BIOMASS LOSS (KG/M2)
!    *PBGAIN*       GRID AVERAGED BIOMASS GAIN (KG/M2)
!    *PBIOMSTR*     GRID AVERAGED ABOVE GROUND STRUCTURAL BIOMASS (KG/M2)
!    *PBIOMSTR2*    GRID AVERAGED BELOW GROUND STRUCTURAL BIOMASS (KG/M2)
!    *PDHBIOS*      Diagnostic array for biomass (see module yomcdh) (KG/M2)
!    *PDHVEGS*      Diagnostic array for variables per vegetation type
!                   (see module yomcdh) 
  
!     AUTHOR
!     ------
!     M.H. Voogt (KNMI) "C-Tessel"  09/2005 
!     S. LAFONT (ECMWF) "C-TESSEL"  05/2006
!     F. Vana  05-Mar-2015  Support for single precision
!     M. Kelbling and S. Thober (UFZ) 11/6/2020 implemented spatially distributed parameters and
!                                               use of parameter values defined in namelist
!     I. Ayan-Miguez (BSC) Sep 2023 Added PSSDP2 object for spatially distributed parameters
!-------------------------------------------------------------------------------

USE PARKIND1  ,ONLY : JPIM, JPRB
USE YOMHOOK   ,ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_CST   ,ONLY : TCST
USE YOS_VEG   ,ONLY : TVEG
USE YOS_AGS   ,ONLY : TAGS
USE NITRO_DECLINE_MOD
USE YOMSURF_SSDP_MOD

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KLON
INTEGER(KIND=JPIM),INTENT(IN)    :: KLEVS
INTEGER(KIND=JPIM),INTENT(IN)    :: KSTEP
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSPHY
REAL(KIND=JPRB)   ,INTENT(IN)    :: PLAT(:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCVT(:,:)
INTEGER(KIND=JPIM),INTENT(INOUT) :: KVEG(:,:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSKM1M(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSAM1M(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PANFMVT(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PANDAYVT(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLAIVT(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PRESPBSTR(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PRESPBSTR2(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PBIOMASS_LAST(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PBIOMASSTR_LAST(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PBIOMASSTR2_LAST(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PBLOSSVT(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PBGAINVT(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLAIE1(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PBSTRE1(:,:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PBSTR2E1(:,:)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PLAI(:)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PBIOM(:)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PBLOSS(:)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PBGAIN(:)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PBIOMSTR(:)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PBIOMSTR2(:)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PDHBIOS(:,:,:)
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PDHVEGS(:,:,:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PLAILC(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PLAIHC(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLAIL(:)
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PLAIH(:)

INTEGER(KIND=JPIM),INTENT(IN)    :: ITVL(:)
INTEGER(KIND=JPIM),INTENT(IN)    :: ITVH(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCVL(:)  
REAL(KIND=JPRB)   ,INTENT(IN)    :: PCVH(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PMU0(:)

REAL(KIND=JPRB)   ,INTENT(IN)    :: PSSDP2(:,:)

TYPE(TCST)        ,INTENT(IN)    :: YDCST
TYPE(TVEG)        ,INTENT(IN)    :: YDVEG
TYPE(TAGS)        ,INTENT(IN)    :: YDAGS

!*         0.     LOCAL VARIABLES.
!                 ----- ----------

INTEGER(KIND=JPIM) :: JL, JVT, IVEG, IVEGTIME

REAL(KIND=JPRB)    :: ZBIOMASSVT(KLON,YDVEG%NVTILES),ZLAI_LAST(KLON,YDVEG%NVTILES),&
                    & ZBIOMASSTR_LAST(KLON,YDVEG%NVTILES),ZBIOMASSTR2_LAST(KLON,YDVEG%NVTILES)
REAL(KIND=JPRB)    :: ZLAICVT(KLON,YDVEG%NVTILES)
REAL(KIND=JPRB)    :: ZXM,ZBMCOEF
REAL(KIND=JPRB)    :: ZEPSILON
REAL(KIND=JPRB)    :: ZLAIMIN
REAL(KIND=JPRB)    :: ZVBSLAI_NITRO
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SURFVEG_EVOL_MOD_CTL:SURFVEG_EVOL_CTL',0,ZHOOK_HANDLE)
ASSOCIATE(RMC=>YDAGS%RMC, RMCO2=>YDAGS%RMCO2, RPCCO2=>YDAGS%RPCCO2, &
 & RVBSLAI_NITROL2D=>PSSDP2(:,SSDP2D_ID%NRVBSLAI_NITROL2D), RVBSLAI_NITROH2D=>PSSDP2(:,SSDP2D_ID%NRVBSLAI_NITROH2D), &
 & RVLAIMINL2D=>PSSDP2(:,SSDP2D_ID%NRVLAIMINL2D), RVLAIMINH2D=>PSSDP2(:,SSDP2D_ID%NRVLAIMINH2D), &
 & RDAY=>YDCST%RDAY, &
 & NVTILES=>YDVEG%NVTILES, RLAIINT=>YDVEG%RLAIINT, RVCOVH2D=>PSSDP2(:,SSDP2D_ID%NRVCOVH2D), &
 & RVCOVL2D=>PSSDP2(:,SSDP2D_ID%NRVCOVL2D))

ZEPSILON= 10._JPRB**(-MAXEXPONENT(ZEPSILON)/10)

IVEGTIME=NINT(RDAY/PTSPHY) !timestep for the vegetation evolution
KVEG(KIDIA:KFDIA,:)=0_JPIM
KVEG(KIDIA:KFDIA,1)=ITVL(KIDIA:KFDIA) !type low veg
KVEG(KIDIA:KFDIA,2)=ITVH(KIDIA:KFDIA) !type high veg KVEG
! SAME FOR VEG COVER
PCVT(KIDIA:KFDIA,:)=0.0_JPRB
PCVT(KIDIA:KFDIA,1)=PCVL(KIDIA:KFDIA)
PCVT(KIDIA:KFDIA,2)=PCVH(KIDIA:KFDIA)
!PCVT(KIDIA:KFDIA,1)=PCVL(KIDIA:KFDIA)/RVCOV(ITVL(KIDIA:KFDIA))
!PCVT(KIDIA:KFDIA,2)=PCVH(KIDIA:KFDIA)/RVCOV(ITVH(KIDIA:KFDIA))
! GOES FROM  LAIL,LAIH to  LAIVT 
PLAIVT(KIDIA:KFDIA,:)=0.0_JPRB
PLAIVT(KIDIA:KFDIA,1)=PLAIL(KIDIA:KFDIA)
PLAIVT(KIDIA:KFDIA,2)=PLAIH(KIDIA:KFDIA)

ZLAICVT(KIDIA:KFDIA,:)=0.0_JPRB
ZLAICVT(KIDIA:KFDIA,1)=PLAILC(KIDIA:KFDIA)
ZLAICVT(KIDIA:KFDIA,2)=PLAIHC(KIDIA:KFDIA)

DO JL=KIDIA,KFDIA
  DO JVT=1,NVTILES
     IF (JVT .EQ. 1_JPIM) THEN
        ! low vegetation
        ZVBSLAI_NITRO = RVBSLAI_NITROL2D(JL)
     ELSE IF (JVT .EQ. 2_JPIM) THEN
        ! low vegetation
        ZVBSLAI_NITRO = RVBSLAI_NITROH2D(JL)
     ELSE
        STOP 'Wrong number of vegitation types in nitro_decline_mod'
     END IF
     ZBIOMASSVT(JL,JVT)=PLAIVT(JL,JVT)*ZVBSLAI_NITRO
     ZBIOMASSVT(JL,JVT)= ZBIOMASSVT(JL,JVT)*(RLAIINT)+ZLAICVT(JL,JVT)*ZVBSLAI_NITRO*(1._JPRB-(RLAIINT))
     ZLAI_LAST(JL,JVT)=PLAIVT(JL,JVT)
     ZBIOMASSTR_LAST(JL,JVT)=PBIOMASSTR_LAST(JL,JVT)
     ZBIOMASSTR2_LAST(JL,JVT)=PBIOMASSTR2_LAST(JL,JVT)
  ENDDO
ENDDO

IF (KSTEP  ==  0) THEN
  DO JL=KIDIA,KFDIA 
    DO JVT=1,NVTILES
      PRESPBSTR(JL,JVT)=0.0_JPRB
      PRESPBSTR2(JL,JVT)=0.0_JPRB
      PBIOMASS_LAST(JL,JVT)=ZBIOMASSVT(JL,JVT)
      PBLOSSVT(JL,JVT)=0.0_JPRB
      PBGAINVT(JL,JVT)=0.0_JPRB
    ENDDO
  ENDDO

ENDIF

IF (RLAIINT > 0._JPRB ) THEN 

! vegetation updated before beginning of new day (in agreement with updates in
! updtim1s at the beginning of the new day) 
    ZBMCOEF=RMC/(RMCO2*RPCCO2)


    DO JVT=1,NVTILES

      CALL NITRO_DECLINE(KIDIA,KFDIA,JVT,KVEG(:,JVT),KLON,KSTEP,PTSPHY, &
         & PCVT(:,JVT), &
         & PTSKM1M,PTSAM1M(:,2),PLAT,PLAIVT(:,JVT),&
         & PSSDP2,&
         & YDCST,YDAGS,&
         & PANFMVT(:,JVT),PRESPBSTR(:,JVT),PRESPBSTR2(:,JVT), &
	 & PBIOMASS_LAST(:,JVT),PBIOMASSTR_LAST(:,JVT),PBIOMASSTR2_LAST(:,JVT), &
	 & ZBIOMASSVT(:,JVT),PBLOSSVT(:,JVT))

      DO JL=KIDIA,KFDIA  

      IF (MOD(KSTEP,IVEGTIME) == 0 ) THEN
         IVEG=KVEG(JL,JVT)
         IF (PCVT(JL,JVT) > 0._JPRB .AND. IVEG >0) THEN
           
           ! computation of LAIGAIN are included in surfveg_evol_ctl

            IF (JVT .EQ. 1_JPIM) THEN
               ! low vegetation
               ZLAIMIN = RVLAIMINL2D(JL)
               ZVBSLAI_NITRO = RVBSLAI_NITROL2D(JL)
            ELSE IF (JVT .EQ. 2_JPIM) THEN
               ! low vegetation
               ZLAIMIN = RVLAIMINH2D(JL)
               ZVBSLAI_NITRO = RVBSLAI_NITROH2D(JL)
            ELSE
               STOP 'Wrong number of vegitation types in nitro_decline_mod'
            END IF
            ! change biomass in time due to assimilation of CO2:
            ZXM=MAX((ZLAIMIN*ZVBSLAI_NITRO-ZBIOMASSVT(JL,JVT)),PANDAYVT(JL,JVT)*ZBMCOEF)
            ! protection
            ZXM=MAX(ZXM,0.0_JPRB)  

            ZBIOMASSVT(JL,JVT)=ZBIOMASSVT(JL,JVT)+ZXM
!           ZBIOMASSVT(JL,JVT)= ZBIOMASSVT(JL,JVT)*RLAIINT+ZLAICVT(JL,JVT)*ZVBSLAI_NITRO*(1._JPRB-RLAIINT)
            ! change in LAI in time due to biomass changes:
            PLAIVT(JL,JVT)=ZBIOMASSVT(JL,JVT)/ZVBSLAI_NITRO


            ! reset to zero the daily net assimilation for next day:
            PANDAYVT(JL,JVT)=0._JPRB

            ! gain output 
            PBGAINVT(JL,JVT)=ZXM

        ELSE 
          PANDAYVT(JL,JVT)=0._JPRB
          PLAIVT(JL,JVT)=0._JPRB
        ENDIF
      ENDIF
      ENDDO  

    ENDDO !veg tile loop


!pb with very low value of PBIOMASSTR_LAST
!where (PBIOMASSTR_LAST<1e-80_JPRB) PBIOMASSTR_LAST=0._JPRB
where (PBIOMASSTR_LAST<ZEPSILON ) PBIOMASSTR_LAST=0._JPRB

ENDIF



! DDH diagnostics
IF (SIZE(PDHVEGS) > 0) THEN
  PDHVEGS(KIDIA:KFDIA,1:NVTILES,2)=PLAIVT(KIDIA:KFDIA,1:NVTILES)
ENDIF
IF (SIZE(PDHBIOS) > 0) THEN
  PDHBIOS(KIDIA:KFDIA,1:NVTILES,1)=ZBIOMASSVT(KIDIA:KFDIA,1:NVTILES)
  PDHBIOS(KIDIA:KFDIA,1:NVTILES,2)=PBLOSSVT(KIDIA:KFDIA,1:NVTILES)
  PDHBIOS(KIDIA:KFDIA,1:NVTILES,3)=PBGAINVT(KIDIA:KFDIA,1:NVTILES)
  PDHBIOS(KIDIA:KFDIA,1:NVTILES,4)=PBIOMASSTR_LAST(KIDIA:KFDIA,1:NVTILES)  
  PDHBIOS(KIDIA:KFDIA,1:NVTILES,5)=PBIOMASSTR2_LAST(KIDIA:KFDIA,1:NVTILES)  
ENDIF

! Grid averaged values
PLAI(KIDIA:KFDIA)=0._JPRB
PBIOM(KIDIA:KFDIA)=0._JPRB  
PBLOSS(KIDIA:KFDIA)=0._JPRB
PBGAIN(KIDIA:KFDIA)=0._JPRB
PBIOMSTR(KIDIA:KFDIA)=0._JPRB
PBIOMSTR2(KIDIA:KFDIA)=0._JPRB
DO JVT=1,NVTILES  
  DO JL=KIDIA,KFDIA
    PLAI(JL)=PLAI(JL)+PCVT(JL,JVT)*PLAIVT(JL,JVT)
    PBIOM(JL)=PBIOM(JL)+PCVT(JL,JVT)*ZBIOMASSVT(JL,JVT) 
    PBLOSS(JL)=PBLOSS(JL)+PCVT(JL,JVT)*PBLOSSVT(JL,JVT) 
    PBGAIN(JL)=PBGAIN(JL)+PCVT(JL,JVT)*PBGAINVT(JL,JVT) 
    PBIOMSTR(JL)=PBIOMSTR(JL)+PCVT(JL,JVT)*PBIOMASSTR_LAST(JL,JVT)    
    PBIOMSTR2(JL)=PBIOMSTR2(JL)+PCVT(JL,JVT)*PBIOMASSTR2_LAST(JL,JVT)    
  ENDDO
ENDDO

!write(6,*) 'vegetation_evol bloss',PBLOSS(1)


!PLAIL(KIDIA:KFDIA)=PLAIVT(KIDIA:KFDIA,1)
!PLAIH(KIDIA:KFDIA)=PLAIVT(KIDIA:KFDIA,2)



! tendency of prognostic variables
DO JL=KIDIA,KFDIA
  DO JVT=1,NVTILES 

    PLAIE1(JL,JVT)=PLAIE1(JL,JVT)+(PLAIVT(JL,JVT)-ZLAI_LAST(JL,JVT))/PTSPHY  

    PBSTRE1(JL,JVT)=PBSTRE1(JL,JVT)+(PBIOMASSTR_LAST(JL,JVT)- &
     & ZBIOMASSTR_LAST(JL,JVT))/PTSPHY
    PBSTR2E1(JL,JVT)=PBSTR2E1(JL,JVT)+(PBIOMASSTR2_LAST(JL,JVT)- &
     & ZBIOMASSTR2_LAST(JL,JVT))/PTSPHY
  ENDDO
ENDDO
END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('SURFVEG_EVOL_MOD_CTL:SURFVEG_EVOL_CTL',1,ZHOOK_HANDLE)

END SUBROUTINE SRFVEGEVOL
END MODULE SRFVEGEVOL_MOD
